Giga Games 1/Giga Games.iso/net/usenet/volume5/pi_perl

Internet Message Format | 1988-08-30 | 9.9 KB

Path: uunet!tektronix!tekgen!tekred!games
From: [e-mail]
Newsgroups: comp.sources.games
Subject: v05i058:  pi-perl - compute the value of pi using perl (and C)
Message-ID: <[e-mail]>
Date: 31 Aug 88 17:24:10 GMT
Sender: [e-mail]
Lines: 418
Approved: [e-mail]

Submitted by: uunet.uu.net!unido!nixpbe!kebsch
Comp.sources.games: Volume 5, Issue 58
Archive-name: pi-perl


    [Is this a game or not? I'm not sure. In any event, it was sent
     to me, so here it is.  -br]

[[The following two programs are for your entertainment. pi.pl is a perl script
an pi.c is the equivalent written in C. Now you may compute pi (3,14....)
about thousands of digits. Yes, the programs are just for fun, especially
the pi.pl (Perl version)! :-)]]

#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create:
#    pi.pl
#    pi.c
# This archive created: Sat Aug  6 12:33:15 1988
export PATH; PATH=/bin:/usr/bin:$PATH
if test -f 'pi.pl'
then
    echo shar: "will not over-write existing file 'pi.pl'"
else
sed 's/^X//' << \SHAR_EOF > 'pi.pl'
Xeval "exec nice -19 perl $0 $*"
X    if $running_under_some_shell;
X# ---------------------------------------------------------------------------
X# pi.perl  computes pi (3.14...) about 5120 Digits
X#
X# W. Kebsch, July-1988  {uunet!mcvax}!unido!nixpbe!kebsch
X
X$my_name = `basename $0`; chop($my_name);
X$version = $my_name . "-1.2";
X
X# some working parameter
X
X$smax =  5120;          # max digits
X$lmax =     4;          # digits per one array element
X$hmax = 10000;          # one array element contains: 0..9999
X$smin = $lmax;          # min digits
X$mag  =     7;          # magic number
X
X# subroutines
X
Xsub mul_tm              # multiply the tm array with a long value
X{
X    $cb = pop(@_);      # elements(array)
X    $x  = pop(@_);      # value
X
X    $c = 0;
X    for($i = 1; $i <= $cb; $i++)
X    {
X    $z      = $tm[$i] * $x + $c;
X    $c      = int($z / $hmax);
X    $tm[$i] = $z - $c * $hmax;
X    }
X}
X
Xsub mul_pm              # multiply the pm array with a long value
X{
X    $cb = pop(@_);      # elements(array)
X    $x  = pop(@_);      # value
X
X    $c = 0;
X    for($i = 1; $i <= $cb; $i++)
X    {
X    $z      = $pm[$i] * $x + $c;
X    $c      = int($z / $hmax);
X    $pm[$i] = $z - $c * $hmax;
X    }
X}
X
Xsub divide              # divide the tm array by a long value
X{
X    $cb = pop(@_);      # elements(array)
X    $x  = pop(@_);      # value
X
X    $c = 0;
X    for($i = $cb; $i >= 1; $i--)
X    {
X    $z      = $tm[$i] + $c;
X    $q      = int($z / $x);
X    $tm[$i] = $q;
X    $c      = ($z - $q * $x) * $hmax;
X    }
X}
X
Xsub add                 # add tm array to pm array
X{
X    $cb = pop(@_);      # elements(array)
X
X    $c = 0;
X    for($i = 1; $i <= $cb; $i++)
X    {
X    $z = $pm[$i] + $tm[$i] + $c;
X    if($z >= $hmax)
X    {
X        $pm[$i] = $z - $hmax;
X        $c      = 1;
X    }
X    else
X    {
X        $pm[$i] = $z;
X        $c      = 0;
X    }
X    }
X}
X
X$m0 = 0; $m1 = 0; $m2 = 0;
X
Xsub check_xb            # reduce current no. of elements (speed up!)
X{
X    $cb = pop(@_);      # current no. of elements
X
X    if(($pm[$cb] == $m0) && ($pm[$cb - 1] == $m1) && ($pm[$cb - 2] == $m2))
X    {
X    $cb--;
X    }
X    $m0 = $pm[$cb];
X    $m1 = $pm[$cb - 1];
X    $m2 = $pm[$cb - 2];
X    $cb;
X}
X
Xsub display             # show the result
X{
X    $cb = pop(@_);      # elements(array);
X
X    printf("\n%3d.", $pm[$cb]);
X    $j = $mag - $lmax;
X    for($i = $cb - 1; $i >= $j; $i--)
X    {
X    printf(" %04d", $pm[$i]);
X    }
X    print "\n";
X}
X
Xsub the_job             # let's do the job
X{
X    $s = pop(@_);       # no. of digits
X
X    $s  = int(($s + $lmax - 1) / $lmax) * $lmax;
X    $b  = int($s / $lmax) + $mag - $lmax;
X    $xb = $b;
X    $t  = int($s * 5 / 3);
X
X    for($i = 1; $i <= $b; $i++)         # init arrays
X    {
X    $pm[$i] = 0;
X    $tm[$i] = 0;
X    }
X    $pm[$b - 1] = $hmax / 2;
X    $tm[$b - 1] = $hmax / 2;
X
X    printf("digits:%5d, terms:%5d, elements:%5d\n", $s, $t, $b);
X    for($n = 1; $n <= $t; $n++)
X    {
X    printf("\r\t\t\t  term:%5d", $n);
X    if($n < 200)
X    {
X        do mul_tm((4 * ($n * $n - $n) + 1), $xb);
X    }
X    else
X    {
X        do mul_tm((2 * $n - 1), $xb);
X        do mul_tm((2 * $n - 1), $xb);
X    }
X    if($n < 100)
X    {
X        do divide(($n * (16 * $n + 8)), $xb);
X    }
X    else
X    {
X        do divide((8 * $n), $xb);
X        do divide((2 * $n + 1), $xb);
X    }
X    do add($xb);
X    if($xb > $mag)
X    {
X        $xb = do check_xb($xb);
X    }
X    }
X    do mul_pm(6, $b);
X    do display($b);
X    ($user,$sys,$cuser,$csys) = times;
X    printf("\n[u=%g  s=%g  cu=%g  cs=%g]\n",$user, $sys, $cuser, $csys);
X}
X
X# main block ----------------------------------------------------------------
X
X$no_of_args = $#ARGV + 1;
Xprint("$version, ");
Xdie("usage: $my_name <no. of digits>") unless($no_of_args == 1);
X$digits = int($ARGV[0]);
Xdie("no. of digits out of range [$smin\..$smax]")
X                unless(($digits >= $smin) && ($digits <= $smax));
Xdo the_job($digits);
Xexit 0;
X
X# That's all ----------------------------------------------------------------
SHAR_EOF
if test 3802 -ne "`wc -c < 'pi.pl'`"
then
    echo shar: "error transmitting 'pi.pl'" '(should have been 3802 characters)'
fi
fi
if test -f 'pi.c'
then
    echo shar: "will not over-write existing file 'pi.c'"
else
sed 's/^X//' << \SHAR_EOF > 'pi.c'
X/*   ------------------------------------------------------------------
X *   Function: Compute a long PI-constant     (4..15,000 see    "smax")
X *
X *   Made by Waldemar Kebsch, Salzkotten, Germany - December 1985
X *                              UUCP: {uunet,mcvax}!unido!nixpbe!kebsch
X *
X *   Usage: pi [number]          (number of digits    on the right hand)
X */
X
X#include <stdio.h>;
X
X#define VERSION  ("  V-1.10  ")
X
X#define     hmax      10000l        /* one element contains: 0..9999 */
X#define     lmax           4        /* max digits of an    element         */
X#define     smin        lmax        /* min digits on the right hand  */
X#define     smax        15000        /* max digits ..             */
X#define     ml   (1+smax/lmax+7-lmax)  /* max. necessary matrix length  */
X
X
Xmain(argc, argv)  /* Parameter:    Number of digits */
Xint argc;
Xchar **argv;
X{          /* 1.    Check Parameter    */
X    int    s;
X
X    if(argc < 2)
X    {
X    printf("Missing argument.\n");
X    exit(1);
X    }
X    if(argc > 2)
X    {
X    printf("Too many arguments.\n");
X    exit(1);
X    }
X    if(sscanf(*++argv, "%d", &s) != 1)
X    {
X    printf("Syntax error.\n");
X    exit(1);
X    }
X    if(s < smin    || s > smax)
X    {
X    printf("Sorry, out of range.\n");
X    exit(1);
X    }
X
X    doit(s);      /* 2.    Let's do it */
X}
X
X
Xdoit(s)       /* Work */
Xint s;
X{
X    char *malloc();
X    long  *pm;           /* primarily matrix */
X    long  *tm;           /* secondary matrix */
X    int      b, xb;       /* number of necessary blocks (elem.) in    the matrix */
X    long  n, t;           /* n=current term, t=max. terms */
X
X    s =    ((s + lmax - 1)    / lmax)    * lmax;        /* digits */
X    b =    xb = s / lmax +    7 - lmax;        /* blocks */
X    t =    s + 2l * s / 3l;            /* terms  */
X    pm = (long*)malloc(b * sizeof(long));   /* array pm   */
X    tm = (long*)malloc(b * sizeof(long));   /* array tm   */
X
X    initial(pm,    tm, b);
X
X    printf("\n");
X    printf VERSION;
X    printf("digits: %5d,  terms: %5d, ", s, t);
X    printf("blocks: %5d\n", b);
X    nice(19);
X
X    for(n = 1l;    n <= t;    n++)   /* time eating loop */
X    {
X    printf("\rterm: %5d ", n); fflush(stdout);
X    if(n < 232l)
X        multiply(tm, (4l * (n * n -    n) + 1l), xb);
X    else
X    {
X        multiply(tm, (2l * n - 1l),    xb );
X        multiply(tm, (2l * n - 1l),    xb );
X    }
X    if(n < 115l)
X        divide(tm, (n * (16l * n + 8l)), xb);
X    else
X    {
X        divide(tm, (8l * n), xb );
X        divide(tm, (2l * n + 1l), xb );
X    }
X
X    add(pm,    tm, xb);
X    xb = checkxb(pm, xb);
X    }
X
X    multiply(pm, 6l, b);
X    display(pm,    b);
X}
X
Xinitial(cm1, cm2, cb)       /* cm1 = cm2    = 0.5 */
Xregister long *cm1, *cm2;  /* current matrix */
Xint cb;
X{
X    long *cmm;
X
X    for(cmm = cm1 + cb;    cm1 < cmm; )
X    *++cm1 = *++cm2    = 0;
X
X    *--cm1 = *--cm2 = hmax / 2l;
X}
X
Xcheckxb(cm, cb)       /* Reduce current block number (speed up) */
Xregister long *cm;
Xint cb;
X{
X    static long    mm = 0,    mn = 0,    nn = 0;
X    register long *cmm;
X
X    cmm    = cm;
X    cm += cb;
X    if(*cm == mm && *--cm == mn    && *--cm == nn)
X    cb--;
X
X    cmm    += cb;
X    mm = *cmm;
X    mn = *--cmm;
X    nn = *--cmm;
X    return(cb);
X}
X
Xmultiply(cm, x,    cb)  /*    cm *= x    */
Xregister long *cm;   /*    current    matrix     */
Xlong x;        /* multiply    with ..    */
Xint cb;        /* current blocks    */
X{
X    long *cmm;
X    long c, z;      /*      carry    */
X
X    for(c = 0, cmm = cm    + cb; cm < cmm;    )
X    {
X    z = *++cm * x +    c;
X    c = z /    hmax;
X    *cm = z    - c * hmax;      /* the same as "*cm = z % hmax;" */
X    }                  /* but faster! */
X}
X
Xdivide(cm, x, cb)  /* cm /= x */
Xregister long *cm;   /*    current    matrix */
Xlong x;        /* divide by ..   */
Xint cb;        /* current blocks */
X{
X    long *cmm;
X    long c, z, q;      /* c=carry, z, q=temporary variable */
X
X    for(c = 0, cmm = cm, cm += cb; cm >    cmm; )
X    {
X    z = *cm    + c;
X    *cm-- =    q = z /    x;
X    c = (z - q * x)    * hmax;      /* the same as  "c = (z % x) * hmax;"    */
X    }                  /* but faster! */
X}
X
Xadd(cm1, cm2, cb)    /* cm1 += cm2 */
Xregister long *cm1, *cm2;    /*    current    matrix */
Xint cb;            /* current blocks */
X{
X    long *cmm;
X    long c;  /*    carry */
X
X    for(c = 0, cmm = cm1 + cb; cm1 < cmm; )
X    {
X    if((*++cm1 += (*++cm2 +    c)) >= hmax)
X    {
X        *cm1 -= hmax;
X        c =    1;
X    }
X    else
X        c =    0;
X    }
X}
X
Xdisplay(cm, cb)
Xregister long *cm;  /* current matrix */
Xint cb;
X{
X    int    i;
X
X    printf("\n%1d. \n",    *(cm +=    cb));     /* one    digit on the left hand */
X    for(i = cb-1; i >= (7-lmax); i--)
X    printf(" %04d",    *--cm);          /* print all digits on the rigth hand */
X
X    printf("\n");
X}
X
X/* That's all */
SHAR_EOF
if test 4336 -ne "`wc -c < 'pi.c'`"
then
    echo shar: "error transmitting 'pi.c'" '(should have been 4336 characters)'
fi
fi
#    End of shell archive
exit 0